/*
   FieldVis - A Tool for Visualizing Astrophysical Magnetohydrodynamic Flow
   
   Copyright (C) 2006 - 2008     FieldVis Development Team
  
   This program is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 2 of the License, or
   (at your option) any later version.
  
   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.
  
   You should have received a copy of the GNU General Public License along
   with this program; if not, write to the Free Software Foundation, Inc.,
   51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 
   File     : CubicSplinePolynomial3D.cpp
   Author   : Blayne Field
   
   Purpose  : Implementation file for the CubicSplinePolynomial class
*/

#include "CubicSplinePolynomial3D.h"


CubicSplinePolynomial3D::CubicSplinePolynomial3D() { }

CubicSplinePolynomial3D::CubicSplinePolynomial3D(Point p0, Point p1,
                                                 Vector3D d0, Vector3D d1)
{
   float distance = Point::DistanceBetween(p0, p1);
   d0.Normalize();
   d1.Normalize();

   ComputeCoefficients(&xCubic, &xQuadratic, &xLinear, &xConstant,
                       p0.GetX(), p1.GetX(), distance * d0.GetX(), 
                       distance * d1.GetX());
   ComputeCoefficients(&yCubic, &yQuadratic, &yLinear, &yConstant,
                       p0.GetY(), p1.GetY(), distance * d0.GetY(),
                       distance * d1.GetY());
   ComputeCoefficients(&zCubic, &zQuadratic, &zLinear, &zConstant,
                       p0.GetZ(), p1.GetZ(), distance * d0.GetZ(), 
                       distance * d1.GetZ());
}

Point CubicSplinePolynomial3D::EvaluateAt(float t) const
{
   // evaluate the polynomial using Horner's rule
   float x = xConstant + t*(xLinear + t*(xQuadratic + t*xCubic));
   float y = yConstant + t*(yLinear + t*(yQuadratic + t*yCubic));
   float z = zConstant + t*(zLinear + t*(zQuadratic + t*zCubic));

   // this should never happen during the execution of the program
   if (t < 0 || t > 1)
   {
      std::cerr << "WARNING: " << t << " is an invalid value for spline evaluation"
                << std::endl;
   }
   return Point(x, y, z);
}

// gives an estimate of the curvature of the spline
float CubicSplinePolynomial3D::AverageCurvature() const
{
   float mag, magCubed;
   float curvature = 0;
   Vector3D r1, r2;
   CubicSplinePolynomial3D firstDeriv = Derivative();
   CubicSplinePolynomial3D secondDeriv = firstDeriv.Derivative();

   // at t = 0
   r1 = firstDeriv.EvaluateAt(0.0).PointToVector3D();
   r2 = secondDeriv.EvaluateAt(0.0).PointToVector3D();
   mag = r1.Norm();
   magCubed = mag * mag * mag;
   curvature += (r1.Cross(r2).Norm()) / magCubed;

   // at t = 0.5
   r1 = firstDeriv.EvaluateAt(0.5).PointToVector3D();
   r2 = secondDeriv.EvaluateAt(0.5).PointToVector3D();
   mag = r1.Norm();
   magCubed = mag * mag * mag;
   curvature += (r1.Cross(r2).Norm()) / magCubed;

   // at t = 1.0
   r1 = firstDeriv.EvaluateAt(1.0).PointToVector3D();
   r2 = secondDeriv.EvaluateAt(1.0).PointToVector3D();
   mag = r1.Norm();
   magCubed = mag * mag * mag;
   curvature += (r1.Cross(r2).Norm()) / magCubed;

   return(curvature / 3);
}

CubicSplinePolynomial3D CubicSplinePolynomial3D::Derivative() const
{
   CubicSplinePolynomial3D deriv;

   deriv.xCubic = 0;
   deriv.xQuadratic = 3 * xCubic;
   deriv.xLinear = 2 * xQuadratic;
   deriv.xConstant = xLinear;

   deriv.yCubic = 0;
   deriv.yQuadratic = 3 * yCubic;
   deriv.yLinear = 2 * yQuadratic;
   deriv.yConstant = yLinear;

   deriv.zCubic = 0;
   deriv.zQuadratic = 3 * zCubic;
   deriv.zLinear = 2 * zQuadratic;
   deriv.zConstant = zLinear;

   return deriv;
} 

void CubicSplinePolynomial3D::ComputeCoefficients(float *cube, float *quad,
                                                  float *linr, float *cons,
                                                  float x0, float x1,
                                                  float dx0, float dx1)
{
   *cube = (2*x0) - (2*x1) + dx0 + dx1;
   *quad = -(3*x0) + (3*x1) - (2*dx0) - dx1;
   *linr = dx0;
   *cons = x0;
}

